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Abstract. A model of unsteady nitration (seepage) in a porous 
medium with capillary retention is considered. It leads to a free 
boundary problem for a generalized porous medium equation where 
the location of the boundary of the water mound is determined as 
part of the solution. The numerical solution of the free boundary 
problem is shown to possess self-similar intermediate asymptotics. 
On the other hand, the asymptotic solution can be obtained from 
a non-linear boundary value problem. Numerical solution of the 
resulting eigenvalue problem agrees with the solution of the partial 
differential equation for intermediate times. In the second part of 
the work, we consider the problem of control of the water mound 
extension by a forced drainage. 



1. Introduction. 

In the present work two problems from the theory of filtration through 
a horizontal porous stratum are considered. First we study a short, but 
intense, flooding followed by natural outflow through the vertical face of 
an aquifer. Further, we consider the possibility to control the spreading 
of the water mound by use of forced drainage at the boundary 

An important practical example of such a problem is groundwater 
mound formation and extension following a flood, after a breakthrough 
of a dam, when water (possibly contaminated) enters and then slowly 
extends into a river bank. 

Consider an aquifer that consists of a long porous stratum with an 
impermeable bed at the bottom and a permeable vertical face on one 
side (Fig. [[]). The space coordinate x is directed along the horizontal 
axis with x = at the vertical face. A water reservoir is located in 
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the region x < 0. We assume that the flow is homogeneous in the y- 
direction. The height of the resulting mound is denoted by z = h(x, t). 
The initial level of water in the stratum is assumed to be negligible. 

The problem is formulated as follows. At some time t = — t < 0, 
the water level at the wall begins to rise rapidly, and water enters 
the porous medium. By time t — 0, the water level at the vertical face 
returns to the initial one. We assume that the distribution at time t — 
is given by h(x, 0) = ho(x) and is concentrated over a finite region [0, d] 
(compactly supported). We also assume that ho(x) is concave down and 
hl(x) is gently sloping. 

In problem 1, water naturally seeps through the boundary back into 
the reservoir, giving the boundary condition h(0,t) = 0. Inclusion of 
the effects of capillary retention into the model distinguishes our case 
from the well know dipole-type problem. The numerical and asymp- 
totic solutions for the source-type boundary conditions were obtained 
in [6]. Most recently, the dipole-type problem with capillary retention 
was studied numericaly and analytically, using Lie-group techniques, 
by B. Wagner in [7]. 

Analysis and numeric computations show that in the case of natu- 
ral outflow, the water mound is not extinguished in finite time. The 
outflow rate cannot be further increased by lowering the level at the 
boundary. In problem 2, in order to control the spreading of the water 
mound, forced drainage is introduced. The problem formulation was 
proposed in [3], where a complete mathematical derivation and rigorous 
analysis can be found. The forced drainage can be implemented, for 
instance, by drilling a number of holes or horizontal wells near the im- 
permeable bottom. In this way, an additional discharge rate is created, 
and the fluid level becomes zero on some interval [0, xi(t)]. 

These are certainly highly idealized problems, but their solutions al- 
low one to extract the qualitative properties and to check the numerical 
methods in solving more realistic problems. 

2. Porous medium equation with capillary retention. 

In the case of seepage and gently sloping profile h 2 (x,t) and in the 
absence of capillary retention, the model of flow in a porous stratum is 
described by the Boussinesq equation ([4] see also [2],[1]): 

(2.1) d t h = nd xx h 2 . 

Here k = kpg/2pm, k is the permeability of the medium, m its poros- 
ity (the fraction of the volume in the stratum which is occupied by 
the pores), p the fluid density, p its dynamic viscosity, and g the ac- 
celeration of gravity. According to the hydrostatic law, water pressure 
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p = pg(h — z). Then, the total head H = p + pgz = pgh is constant 
throughout the height of the mound. Under the assumption of seepage 
and gently sloping profiles h 2 , Darcy law is used to obtain the relation 
for the total flux q = —-VH ■ h. 

Mathematical properties of the Boussinesq equation are well known [5]. 
An essential feature of this equation is the finite speed of disturbance 
propagation given a finite (compactly supported) initial distribution. 
Another important feature of this equation is the existence of special 
self-similar solutions. The graphs of such a solution for any two times to 
and t\ are related via a similarity transformation [1]. The special solu- 
tions, themselves corresponding to certain, sometimes artificial, initial 
and boundary conditions, are important because they provide inter- 
mediate asymptotics for a wide class of initial value problems. For 
these problems, the details of the initial distribution affect the solu- 
tion only in the beginning; after some time, the solution approaches a 
self-similar asymptotics. The Boussinesq equation has been studied ex- 
tensively and a number of self-similar solutions, for different boundary 
conditions, have been constructed ([2], [3]). 

Following [1], [6], the Boussinesq equation can be modified to incor- 
porate the effects of capillary retention into the model. If we exclude 
the possibility of water reentering the region that was filled with water 
at some earlier time and assume that initially the stratum is empty, 
we have the following situation: when water enters a pore, it occupies 
the entire volume, allowed by active porosity; when water leaves the 
pore, a fraction S of the pore volume remains occupied by the trapped 
water. We assume that 5 is constant. Let us denote the initial active 
porosity by m. Then, when water is entering previously unfilled pores, 
the effective porosity is m; when water is leaving previously water-filled 
pores, the effective porosity becomes m(l — 5). Hence, in the presence 
of capillary retention, porosity depends on the sign of d t h(x, t). Notice, 
that permeability can be assumed unaffected, as the effect of capillary 
forces on permeability is significant only for small and/or dead-end 
pores, whose contribution to the total flux, in the first approximation, 
can be neglected. 

The rate of change in the amount of water AV inside a volume 
element (Fig. [D) is equal to: 



(2.2) 
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On the other hand, the rate of change in the volume of water due to 
the flux through the faces of a volume element (x, x + Ax) is equal to 



(2.3) Ax d x (-d x H -h) = Ax ^-d xx (h 2 ). 

fj, 2/j, 

We denote «i = jf^ and k 2 = 2l im{i-8) • Then, using the continuity of 
flux (no sources inside the water mound) and the balance of mass we 
obtain: 

(24) dh = U 1 d xx (h 2 ) iff>o, 

{ ' } at \K 2 d xx {h 2 ) iff<o. 

This is a nonlinear parabolic partial differential equation with discon- 
tinuous coefficients, also known as the generalized porous medium equa- 
tion [2], [6]. 

Continuity of the flux q = —P&h ■ d x h implies that at the mound 
tip x r , where mound height is zero, the flux is also zero. For prob- 
lem 1, these considerations lead to the following initial and boundary 
conditions to supplement equation (|2.4|): 



h(x, 0) = h (x) > (where Hq{x) = for x > d), 
[2.5) h(x,t) = 0, d x h 2 (x, t) = at x = x r , 
h(0,t) = 0. 



The second line in ( |2.5| ) corresponds to the free boundary conditions 
on the right boundary, x r (t), which is unknown a priori. 

It should be noted that for the solution of equation Q2.1D (but not 



for (|2.4j) ) with boundary conditions ( |2.5|) the dipole moment is con- 



stant: 



(2.6) Q = J xh(x,t)dx = C. 



We call equation (|2.4j ) with boundary conditions (|2.5f) a dipole-type 
problem. A similar problem, for source type initial and boundary con- 
ditions was considered in [6], see also [1]. 

For problem 2, the boundary conditions are changed to include the 
forced drainage condition. The discharge rate qo(t), which is a quantity 
that should be specified, determines the boundary condition at the left 



free boundary x\. 



h(x, 0) = ho(x) > (where ho(x) = for x > d), 
[2.7) h(x,t) = 0, d x h (x, t) = at x = x r , 

h(x,t) = 0, d x h 2 (x,t) = ^ - at x = xi . 

run 



The second and third lines in ( |2.7| ) define, respectively, the free bound- 
ary condition on the right boundary and the forced drainage condition 
on the left boundary. Equation (|2.4|) together with boundary condi- 
tions ( |2.7| ) define problem 2. 

3. Dimensional analysis of problem 1. 

The parameters in the problem are h, x, t, Ki, k>2, d - the initial width 
of the water mound, and Q = Q(0) - the initial dipole moment. We 
can take the dimensions as follows: [h] = H, [x] = L, [t] = T. Then 
from equation (|2.4j ) we have [k] = t^. For the remaining parameters 
[d] = L, [Q] = H ■ L 2 . The dimensions for h and x are set to be 
independent. This can be done because the differential equation ( |2.4| ) 
is invariant with respect to the following group of transformations: 

a 2 

(3.1) x' = ax, t = — t, h = 'yh. 
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The invariance insures that we can scale the units of measurement for 
h, while keeping the units for x unchanged. 

The following dimensionless quantities can be obtained from these 
parameters: 

n ' = ww^> = n^f^ = S' md n = hC i ,1/2 - 

it follows that n = F(rii, n 2 , n 3 ). 

Since for large times, t 3> the parameter II2 <C 1, it would seem 
natural to set II = /(IT, II3), as in the case of K\ = K2, and look for a 
solution of the form: 

(3.2) h=(^)V 2 f(z,^), where z 



Kit JK 'fc 2 7 ' (QM) 1/4 ' 

However, this leads to a contradiction when we consider an ordinary 
differential equation obtained from (|2.4j ): 



df ., I if 2/ + £ z < 0, 



(3 ' 3) (2/ + ^"i-4|% H2f+ d {z>0. 



Multiplying both sides by z we obtain an equation in total differentials, 
which is readily solved: 



(3.4) fz* 



-4(^1 -m^i if 2/ + £,<(), 



Observe that near z = z r , where the height of the mound vanishes, the 
first equation holds. At z = z r , h vanishes along with the flux, which 
is proportional to From the first equation at z r , we obtain that 
C\ = 0. Similarly, evaluating the second expression at z — 0, where 
h = 0, we find that C 2 = 0. Next, evaluating the two expressions at 
z±, we obtain: 

(3-5) (f^Wi-/ 2 (^i))(l--)=0, 
az K 2 

and using 2/ + |z = we obtain -5/ 2 (zi)(l - — 0. For k 2 = «i, 
the solution can be found [2] and thus the assumption of complete 
similarity in n 2 is correct. However, in the case of n 2 ^ K\, we have 
f(z\) = 0. This is a contradiction, because the change in sign of 8fe g^ 
should occur inside the mound, where the height is positive. Hence, 
the assumption of complete similarity for k 2 ^ n± does not hold. 

We next solve the problem numerically and study the asymptotic 
behavior of the solution. 



4. Numerical solution of the partial differential 
equation and further analysis for problem 1. 



In order to simplify the numerical solution for equation (|2.4j ) with 
free boundary conditions (|2.5| ), we use a change of variables: £ = 
We set H(^,t) = h(x r £,,t), and equation (|2]4]) is transformed: 

(4.1) 

'^( Kl d^H 2 (^ t) + Ki^ff (1, t)d^H(t t)) if dtH(£, t) > 0, 
^{K 2 d^H\i, t) + Kx&tHQ., t)dsH(£, t)) if d t H(£, t) < 0; 



dtH 



with boundary conditions H(0,t) = H(l,t) = 0. This effectively fixes 
the right boundary at £ = 1. 

The location of the free boundary x r (t) can be obtained in the course 
of the numerical solution in the following way. We assume that the 
solution is nearly stationary near the tip x r and h(x,t) ~ H{x — vt). 
Here, v denotes the instantaneous speed of mound extension, which 
changes slowly as a function of t. Then d t h ~ —vd x h near x r . 



Considering equation (|2.4|) near the boundary x r of the mound, 
where h is small, we have: 

d t h = K X d xx h 2 {x,t) = 2K 1 ((d x h(x 1 t)) 2 + hd xx h(x,t)) « 2K 1 (d x h(x,t)) 2 
so that 

(4.2) v(t) = -2Kid x h(x r ,t), x r (t) = [ v(t)dt + x r (0). 

Jo 

We solve the new boundary value problem numerically by using a 
forward-in-time, centered-in-space finite-difference approximation, where 
u™ is an approximation to the solution of (|4.1|) at the grid point (xi, t n ): 



k, if [{vJUlf - 2(ur 1 ) 2 + {u^f\ > 0, 
\k 2 if Mil? - 2(ur 1 ) 2 + K+i 1 ) 2 ] < 0; 

< +1 = < + ^{<[«-i) 2 -2«) 2 + « +1 ) 2 ] 



«i^«-<-i)K-<-i)}/(^ 



n\2 



x — x r 2k\ 



At un — un-i 
Ax x 7 } 



In the numerical computation we start with an initial distribution of 
the source type, localized near x = (Fig.|[). Before the left free 
boundary reaches the point x — 0, the solution is of the source type 
and we can compare our numerical results to those in [6]. After some 
time t the left free boundary reaches x — 0, where it is thereafter fixed 
(Fig. 1). 

Now we consider the scaled solution: 

(4.3) /({| g. )|{= *. 

max£ H (4 j t) k 2 x r 

We can see in Figs. ^ and [7] that as time t increases the numerical solu- 
tion approaches a self-similar regime, so that the graphs of the scaled 
solution for different times "collapse" into a single curve. Moreover, 
Figs. |3] and § show a power-law dependence on time for both maxg H 
and x r in the self-similar regime. 

In part 3, we have shown that a self-similar solution of the first kind 
does not exist for this problem. To explain what happened we return 
to the dimensional analysis and look now for a generalized self-similar 
solution. 

We have determined that the variables in the problem are related as 
follows: IT = F(Ui, U 2 , U 3 ), where 

Hl = Tn^Wf 112 = 7n4wi' U * = ~> and n = h ^"- 
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Our numerical investigation shows that for large t, as II2 — > 



(4.4) F(n 1; n 2 , n 3 ) = iq^u^, n 



3) 



where 7 and e are constants. In fact, this is the next simplest situation 
after complete self-similarity and it is referred to as self-similarity of 
the second kind in n 2 (see [2]). 
Indeed, from the analysis above: 



III 



n 



7 



xd e (tKiQ) 
d? 



5-1 
4 



2 {tK\QY '' 

X K\ . 



(4.5) h(x,t) = Ar«f( B/ r/;>) . 

x r {t) = Bt 13 , 

where 

A = (K l Q)^ /i d" f , B = d £ {KxQ)^,a = and (3 = 

The parameters a and (3 depend on the ratio They cannot be 
determined on the basis of dimensional analysis alone and have to be 
computed as a part of the solution. We will see that there is, actually, 
only one unknown parameter involved, since the differential equation 
provides an additional relation between a and (5. 



5. Derivation and numerical solution of a nonlinear 

eigenvalue problem. 

The numerical solution of partial differential equation showed that 
there is indeed an intermediate asymptotic solution of the form ( |4.5| ). 
Now, we can obtain such a self-similar solution by transforming the 
problem of solving partial differential equation (|2.4| ) with boundary 
conditions ( |2.5| ) into a nonlinear eigenvalue problem. 

We substitute ( |4.5|) into ( |2.4|) and normalizing / so that ^ = 1 we 
get: 

(5.1) dth 

(5.2) d xx h 2 

(5.3) af(0 + 0f(t)t 



= -t-W A(af + /3fZ), 
BH 2 ? 

= -«*(/ 2 (0) ,, *- a - w , 
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where 

K * f«i if ((1-2/3)/ + /3/'0<0, 
\k 2 if ((1-2/3)/ + /?/'£) >0. 

Since equation ( |5~3| ) cannot depend on time explicitly, —a — 2(3 + 1 = 0. 
Finally, we get an ordinary differential equation: 

-{fin + (1 - 2/3)/) if ((1 - 2/3)/ + (3fO < 0, 
-*W£ + (1 - 2/3)/) if ((1 - 2/3)/ + 0f > 0. 



(5-4) (/ 



The boundaries xo = and x r = -Bt^, in the new space variable, 
correspond to £ = and £ = 1. The boundary condition at £ = 
becomes: 

(5.5) /(0)=0. 
For the right boundary, £ = 1, we have: 

(5.6) /(l) = 0and(/ 2 )'(l)=0. 



From g and ^§ it follows that 2/ti(/'(l)) 2 + /3£/'(l) = and the tip 



conditions become: 

(5.7) /(l) = and f'(l) 



2k, 



The second order ODE (|5.4|) with three boundary conditions (5. i 



and |5.5[) constitutes a non-linear eigenvalue problem, which we now 



have to solve numerically. For each value of we find a value of (3 
such that the boundary conditions are satisfied. 

We use a high order, Taylor-expansion-based method to start the 
integration at £ = 1 followed by a 4th order Runge-Kutta method and 
an iterative procedue to arrive at the value for (3 such that the third 
condition /(0) = is satisfied. For computational convenience, we 
transform the differential equation by changing variables: = / 2 (£)> 
so that does not have a singularity at £ = 0. In this manner, we 
obtain the dependence of (3 on ^ (Fig. ffy. 

6. Comparison of the results for problem 1. 
In logarithmic coordinates, we obtain from ( |4.5| ): 

(6.1) log(u(l/2,f)) = -alog(t) + log(A/(l/2,^)) 

(6.2) log(x r (t)) = (3\og(t) + \og(B) 



i.e., straight lines with slopes —a and (3. 
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From plots in Figs. |] and |], we can observe that after some initial 
time t both graphs approach straight lines. We repeat the calculations 
for a range of values of — . 

Comparison of the results of the numerical solution of the nonlinear 
eigenvalue problem with the results obtained from the numerical solu- 
tion to the partial differential equation (Fig. [5]), shows that the two 
agree with high precision. 

Also, the exact solution for the case K\ = K2 gives the value (3 = .25, 
which coincides with the results of the numerical computations with 
good accuracy. 

7. Numerical solution of the partial differential 
equation for problem 2. 

Although problems 1 and 2 are similar, the numerical treatment of 
problem 2 is more complicated. Time evolution of the left boundary in 
problem 2 makes rescaling, which was used in the numerical solution 



of problem 1, infeasible. Instead, we solve equation fl2.4f) on a grid, 
taking into account that the left and right boundaries may not fall 
onto gridpoints. We determine new positions of the boundaries from 
the numerical solution at each timestep. 



Equation (2A) is discretized using a forward- in-time, centered- in- 
space finite-difference scheme: 

Kl if [{utlf - 2K- 1 ) 2 + [u^f] > 0, 

k 2 if [{utlf - 2« -1 ) 2 + K+i 1 ) 2 ] < °; 

< +1 = < + ^R[«-i) 2 - 2«) 2 + « +1 ) 2 ]}, 

(7-1) < +1 = + ^yA iu K~J u?)2 ~ <i n h 

7/ n+l _ ? .n , 2At ,nf «-i) 2 -«) 2 _ «) 2 1 

Here u\ and u r are the nonzero values of u on the grid, adjacent to the 
left and right boundaries respectively, q is the drainage flux, Axi and 
Ax r are distances from the left and right boundaries to the grid points. 
We treat the values ui and u r separately in order to incorporate the 
boundary conditions and improve precision. 

The location of the left boundary is obtained from the values of u: 



(7.2) x? +1 = x? + A Xl - 

The right boundary location is obtained by extrapolation from the 
values of u n+1 . 
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We check the numerical method for K\ = K2 by comparing the nu- 
merical solution with a known analytic solution. The exact self-similar 
solutions for the problems with forced drainage are given in [3]. We 
choose a value of (3 7^ .25 and then solve an ordinary differential equa- 
tion fl5.3|) with the initial condition (|5.7| ). For /3 < .25 the solution /(£) 
of the ordinary differential equation intersects the x-axis at some point 
£ = A > and at £ = 1. 

From the solution of the ordinary differential equation we obtain a 
self-similar solution: 



of the partial differential equation. The locations of the free boundaries 
are given by x r (t) = At 13 , Xi(t) = XAt 13 , and the drainage flux is given 
by q{t) = mBAt~ 2+3l3 (f 2 )'(X) (see [3]). We use the self-similar solution 
at some time to as an initial value for the numerical solver, set drainage 
flux on the left boundary to be q(t), and compute the solutions until 
time t\. As in the analysis in section the graphs of x r (t), Xi(t) 
should be straight lines in logarithmic coordinates, and the graphs of 
the scaled solution for different times t should collapse into one curve. 
That's what we observe in Figs. |9] and |K]. 

Now, we try to model the conditions of a flood followed by forced 
drainage, as described in the introduction. We begin by computing 
the solution to problem 1 until some time to, which corresponds to the 
flood followed by natural drainage through the boundary of the aquifer. 
After to, we set a constant drainage flux q(t) = q at the left boundary. 
In particular, we set q to equal twice the natural drainage flux at time 
to. As we see in Fig. [T2|, the water mound, that has appeared after the 
flood, is completely extinguished in finite time. 



1. The numerical simulations of two problems involving drainage and 
capillary retention of the fluid a in porous medium were presented. 
It was shown that the problem with dipole type initial and bound- 
ary conditions has a self-similar intermediate asymptotics in the 
case of a porous medium with capillary retention. 

2. A problem of control of the water mound extension by forced 
drainage was considered. The possibility of extinguishing the 
propagating water mound by creating a forced drainage flux q(t) 
at the left boundary was confirmed numerically. Using our re- 
sults, it should be possible to derive a cost efficient drilling regime 
and to localize the mound and contain the contamination inside 
a prescribed region. 



(7.3) 




8. Conclusion. 
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It would be interesting to extend the numerical investigation above to 
the case of a fissurized porous medium. 
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Figure 1 . Groundwater dome extension in a porous medium. 
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Figure 3. Example of the time evolution for the posi- 
tion x r of the free boundary. Time interval t = [0,450]. 
For this plot ^ = .5. 
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logt 

Figure 4. Plot of log(max I u(x, t)) vs. log(£) with ^ = 
.5. The straight line shows a linear fit to the straight part 
of the graph. 
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Figure 5. Plot of the dependence of (3 on ^. The solid 
line is obtained from graph |] for different values of 
by applying least-squares fit to the straight part of the 
graph and then using equation (|6.1|) . The dashed line 
is obtained from graphs ||, by applying the procedure 
above and then (|6.2|) . The dotted line is obtained from 
the solution of the nonlinear eigenvalue problem given 
by equation (|5.4|) with initial conditions |x7| for different 
values of — . 

K2 
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Figure 7. Plot of the scaled solution (fOl ) of the partial 
differential equation for time interval t = [18, 450] . For 
this plot ^ = .5. 



19 



0.035r 



0.03- 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



t 

Figure 8. Example of convergence of the "shooting" 
method for solving the nonlinear eigenvalue problem 
given by 0, 
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Figure 9. Positions of the boundaries and max^ h(x, t) 
vs t in logarithmic coordinates for the self-similar regime 
in the forced drainage problem. f3 — .2 
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Figure 10. "Collapse" of the graphs of the scaled so- 
lution in self-similar coordinates for problem 2. Time 
t = [1,10], /3 = .2 
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FIGURE 12. Extinguishing of the water mound with 
forced drainage. 



24 



